import numpy as np
import ipywidgets as ipw

from bokeh.io import push_notebook, output_notebook, show
from bokeh.layouts import row, column
from bokeh.plotting import figure
from bokeh.models import ColumnDataSource
output_notebook()

from ipywidgets import interact

from collections import OrderedDict
old_settings = np.seterr(over = 'ignore') #Ignore warnings about overflow data points

# Set up constants and parameters in SI units
m = 9.109e-31; #Define electron mass
hbar = 6.626e-34/2.0/np.pi; #Reduced planck constant
size = 100e-9; #Define the width of the potential well in real space
Emax = 2.0; # Pick the highest energy (eV) the electron could have by distribution
k0 = np.pi/size; #Lowest allowed k in reciprocal space

# Set up the variable space
N = int(np.sqrt((Emax*1.602e-19*2*m)/hbar**2)/k0); # calculate the number of allowed states by the energy constraint
k = np.linspace(-int(np.sqrt((Emax*1.602e-19*2*m)/hbar**2))/1e9, int(np.sqrt((Emax*1.602e-19*2*m)/hbar**2))/1e9, 2*N); # Double the variable space by symmetry
En = 6.242e+18*(hbar*k*1e9)**2/(2*m);  #Energy for a free electron in eV

# Plot the dispersion relation of the free electrons
freeE = figure(height=350, width=300, title="Dispersion relationship of a free electron",
              tools="pan,reset,save,wheel_zoom", x_range = [np.amin(k), np.amax(k)],
              y_range = [0,1.1*np.amax(En)]);
freeE.xaxis[0].axis_label='k (nm-1)';
freeE.yaxis[0].axis_label='Energy (eV)';
source = ColumnDataSource(data = {'xVal':k, 'yVal':En});
freeE.scatter('xVal', 'yVal', source = source, line_width=1, line_alpha=0.2);

# Plot the histogram of states per energy interval
bWidth = 0.02; # Define the bin width of the histogram as 0.02 eV
histoN = (np.amax(En)/bWidth).astype(int); # calculate the number of bars in the histogram
kCount = np.empty(histoN, dtype = np.int8); # Declare an array to store histogram information
Ens = np.arange(histoN) * bWidth;
zero = np.zeros(histoN);
for idx in range(histoN):
    kCount[idx] = len(En[ (En > idx*bWidth) & (En <= (idx+1)*bWidth) ]);
histo = figure(height=350, width=300, title="allowed states per Energy interval",
              tools="pan,reset,save,wheel_zoom", y_range = [0, 1.1*np.amax(En)],
              x_range = [0, 1.1*np.amax(kCount)]);
histo.xaxis[0].axis_label='Number of states';
histo.yaxis[0].axis_label='Energy (eV)';
source1 = ColumnDataSource(data = {'bottom': Ens, 'top': Ens + bWidth, 'right': kCount, 'left': zero});
histo.quad(bottom = 'bottom', top = 'top', right = 'right', left = 'left', source = source1,
           fill_color="navy", line_color="white", alpha=1);

#calculate Density of states
gE = 1/(np.pi*hbar)*np.sqrt(m/(2*1.602e-19*En));
gEplot = figure(height=350, width=300, title="Density of States (1D)",
              tools="pan,reset,save,wheel_zoom", x_range=[0.5*np.amin(gE), 5*np.amin(gE)], y_range=[0, 1.1*np.amax(En)])
gEplot.xaxis[0].axis_label='Density of States';
gEplot.xaxis[0].formatter.precision = 2;
gEplot.xaxis[0].ticker.desired_num_ticks = 5;
gEplot.yaxis[0].axis_label='Energy (eV)';
source2 = ColumnDataSource(data = {'xVal': gE, 'yVal': En});
gEplot.scatter('xVal', 'yVal', source = source2, line_width=1, line_alpha=0.2);

# Set up callbacks to live update the plots
def update_data(size, Ntot):
    # Generate the new curve
    k0 = np.pi/(size);
    Emax = 2*(hbar*k0*Ntot)**2/m/1.602e-19;
    N = Ntot; # calculate the number of allowed states by the energy constraint
    k = np.linspace(-N*k0, N*k0, 2*N);
    En = 6.242e+18*(hbar*k*1e9)**2/(2*m);
    
    update = {'xVal': k, 'yVal': En};  #Energy for a free electron in eV
    source.data = update;

    histoN = (np.amax(En)/bWidth).astype(int); # calculate the number of bars in the histogram
    kCount = np.empty(histoN, dtype = np.int8); # Declare an array to store histogram information
    Ens = np.arange(histoN) * bWidth;
    zero = np.zeros(histoN)
    for idx in range(histoN):
        kCount[idx] = len(En[ (En > idx*bWidth) & (En <= (idx+1)*bWidth) ]);  
    update1 = {'bottom': Ens, 'top': Ens + bWidth, 'right': kCount, 'left': zero};
    source1.data = update1;

    gE = 1/(np.pi*hbar)*np.sqrt(m/(2*1.602e-19*En))
    update2 = {'xVal': gE, 'yVal': En}
    source2.data = update2;
    push_notebook(handle = handel);

handel = show(row(freeE,histo,gEplot), notebook_handle = True);

interact(update_data, size = ipw.FloatSlider(min = 2, max = 202, step = 10, value = 50, description = 'Size (L, nm)'), Ntot = ipw.IntSlider(min = 4, max = 604, step = 10, value = 200, description = 'No. Electrons'));

interactive(children=(FloatSlider(value=50.0, description='Size (L, nm)', max=202.0, min=2.0, step=10.0), IntS…

This notebook shows density of states in 1D. Note that all graphs are cut off at Fermi Energy Ef, hence if the cut is visible on the graph consider increasing number of electrons in the system.¶